En el flujo con una cavidad rectangular, el fluido se encuentra confinado en una región rectangular de dimensiones $L\times W$ donde el movimiento se debe al desplazamiento transversal de una de las caras. En este sentido, la velocidad de la cara superior es de $U_0$ y el de las demás interfaces es de $0$. De acuerdo con esto, se tienen las siguientes consideraciones:
De aquí, las ecuaciones de continuidad y de Navier-Stokes se reducen a:
\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho \vec{u})=0 \quad &\to\quad \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 \\ \rho\left(\frac{\partial \vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}\right)=\rho\vec{f}-\nabla p+\mu\nabla^2\vec{u}\quad &\to\quad \begin{cases} \rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) \\ \rho \left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial v}+\mu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right) \end{cases} \end{align}Se puede ver entonces que la ecuación resultante es un caso particular del problema de convección-difusión estable:
$$\frac{\partial f}{\partial t} + \vec{u}\cdot\nabla f=k\nabla^2f+\frac{s}{\rho c_p}$$Donde $\frac{\partial f}{\partial t}=0$, $k = \frac{\mu}{\rho}=\nu$, y $\frac{s}{c_p}=-\frac{\partial p}{\partial *}$ para $f=u, *=x$ y $f=v, *=y$, y las condiciones de frontera corresponden a valores conocidos (tipo Dirichlet).
Ahora, para aplicar el método de elementos espectrales, es necesario realizar los siguientes pasos:
Cada uno de los pasos se tratará en detalle a continuación
In [1]:
# Prepara el ambiente de trabajo con las librerías a utilizar
import numpy as np # Ofrece funciones para operaciones básicas con arreglos
import scipy.linalg as lin # Ofrece funciones para operaciones de álgebra lineal
import matplotlib.pyplot as plt # Permite incluir un ambiente de visualización
from mpl_toolkits.mplot3d.axes3d import Axes3D
import timeit # Permite la medición de tiempos para los distintos algoritmos
from sem2D import * # Agrupa las funciones externas en un archivo
# from __future__ import division # Corrige posibles problemas que surgen a la hora de dividir enteros
#%pylab inline
In [2]:
L = 2.4
W = 1.6
mu = 1.
U0 = 1.
Ne = np.array([8,8]) # Indica la cantidad de elementos a lo largo de los ejes x y y
Np = np.array([2,2]) # Indica la cantidad de divisiones sobre cada elemento en xi y en eta
# Genera los puntos correspondientes a los vertices de cada elemento
(px, hx) = np.linspace(-L/2,L/2,Ne[0]+1,retstep=True)
(py, hy) = np.linspace(-W/2,W/2,Ne[1]+1,retstep=True)
In [3]:
(x, y, ne, npe) = genNodes(Ne,Np,px,py,ordn='ij') # Notese que el ordenamiento es ij
# Grafica la malla de puntos
for l in range(ne):
plt.plot([x[l,0],x[l,Np[1]],x[l,-1],x[l,-Np[1]-1],x[l,0]],
[y[l,0],y[l,Np[1]],y[l,-1],y[l,-Np[1]-1],y[l,0]],
'r--')
plt.hold(True)
plt.plot(x,y,'bo',ms=4)
plt.xlabel('coordenada x')
plt.ylabel('coordenada y')
plt.show()
In [4]:
h = min(hx,hy)
(coords, C, gfl, ng) = genConn(Ne,Np,ne,npe,x,y,L,W,U0,h)
In [5]:
(Mat_dif, vec_b) = matDiff(Ne,Np,ne,npe,ng,x,y,C,gfl,mu)
plt.spy(Mat_dif)
plt.show()
Teniendo ya construidos los componentes del sistema a resolver, basta aplicar un esquema de solución razonable para el tipo de matriz obtenida. En este caso se utiliza el solucionador por defecto del módulo linalg. Dado que los valores para el primer y último nodo son conocidos, se omiten los valores correspondientes.
In [6]:
start = timeit.default_timer()
sols = lin.solve(Mat_dif[1:,1:],vec_b[1:])
print timeit.default_timer()-start
# Recupera las velocidades
solx = np.concatenate([[gfl[0,1]], sols[0:ng-1]])
soly = sols[ng-1:2*ng-1]
# Calcula la magnitud de la velocidad en cada nodo
solm = np.abs(solx**2+soly**2)
# Recupera las presiones
solp = sols[2*ng-1:]
In [9]:
# Grafica las fronteras de los elementos
for l in range(ne):
plt.plot([x[l,0],x[l,Np[0]],x[l,-1],x[l,-Np[0]-1],x[l,0]],
[y[l,0],y[l,Np[0]],y[l,-1],y[l,-Np[0]-1],y[l,0]],
'r--')
plt.hold(True)
# Grafica la magnitud del campo de velocidades
for l in range(ne):
Px = (coords[C[l,:],0]).reshape(Np[0]+1,npe/(Np[0]+1))
Py = (coords[C[l,:],1]).reshape(Np[0]+1,npe/(Np[0]+1))
SM = (solm[C[l,:]]).reshape(Np[0]+1,npe/(Np[0]+1))
plt.pcolor(Px,Py,SM,vmin=min(solm), vmax=max(solm))
# Grafica el campo de velocidades
plt.quiver(coords[:,0],coords[:,1],solx,soly,color='w')
plt.title('Campo de Velocidades para el Flujo en la Cavidad')
plt.xlabel('Posicion (coord x)')
plt.ylabel('Posicion (coord y)')
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.show()
In [8]:
# Grafica el campo de presiones
xpres = np.zeros(ne)
ypres = np.zeros(ne)
for l in range(ne):
xpres[l] = (x[l,0]+x[l,Np[0]]+x[l,-1]+x[l,-Np[0]-1])/4
ypres[l] = (y[l,0]+y[l,Np[0]]+y[l,-1]+y[l,-Np[0]-1])/4
xpres = xpres.reshape(Ne[0],Ne[1])
ypres = ypres.reshape(Ne[0],Ne[1])
zpres = solp.reshape(Ne[0],Ne[1])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(xpres,ypres,zpres)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Presion')
#plt.pcolor(xpres,ypres,zpres)
plt.show()
In [ ]: